home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / gems / g_gemsi.lha / GraphicsGems / Sturm / sturm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-10  |  5.1 KB  |  288 lines
  1.  
  2. /*
  3.  * sturm.c
  4.  *
  5.  *    the functions to build and evaluate the Sturm sequence
  6.  */
  7. #include <math.h>
  8. #include <stdio.h>
  9. #include "solve.h"
  10.  
  11. /*
  12.  * modp
  13.  *
  14.  *    calculates the modulus of u(x) / v(x) leaving it in r, it
  15.  *  returns 0 if r(x) is a constant.
  16.  *  note: this function assumes the leading coefficient of v 
  17.  *    is 1 or -1
  18.  */
  19. static int
  20. modp(u, v, r)
  21.     poly    *u, *v, *r;
  22. {
  23.     int        k, j;
  24.     double    *nr, *end, *uc;
  25.  
  26.     nr = r->coef;
  27.     end = &u->coef[u->ord];
  28.  
  29.     uc = u->coef;
  30.     while (uc <= end)
  31.             *nr++ = *uc++;
  32.  
  33.     if (v->coef[v->ord] < 0.0) {
  34.  
  35.  
  36.             for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  37.                 r->coef[k] = -r->coef[k];
  38.  
  39.             for (k = u->ord - v->ord; k >= 0; k--)
  40.                 for (j = v->ord + k - 1; j >= k; j--)
  41.                     r->coef[j] = -r->coef[j] - r->coef[v->ord + k]
  42.                     * v->coef[j - k];
  43.     } else {
  44.             for (k = u->ord - v->ord; k >= 0; k--)
  45.                 for (j = v->ord + k - 1; j >= k; j--)
  46.                 r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  47.     }
  48.  
  49.     k = v->ord - 1;
  50.     while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
  51.         r->coef[k] = 0.0;
  52.         k--;
  53.     }
  54.  
  55.     r->ord = (k < 0) ? 0 : k;
  56.  
  57.     return(r->ord);
  58. }
  59.  
  60. /*
  61.  * buildsturm
  62.  *
  63.  *    build up a sturm sequence for a polynomial in smat, returning
  64.  * the number of polynomials in the sequence
  65.  */
  66. int
  67. buildsturm(ord, sseq)
  68.     int        ord;
  69.     poly    *sseq;
  70. {
  71.     int        i;
  72.     double    f, *fp, *fc;
  73.     poly    *sp;
  74.  
  75.     sseq[0].ord = ord;
  76.     sseq[1].ord = ord - 1;
  77.  
  78.  
  79.     /*
  80.      * calculate the derivative and normalise the leading
  81.      * coefficient.
  82.      */
  83.     f = fabs(sseq[0].coef[ord] * ord);
  84.     fp = sseq[1].coef;
  85.     fc = sseq[0].coef + 1;
  86.     for (i = 1; i <= ord; i++)
  87.             *fp++ = *fc++ * i / f;
  88.  
  89.     /*
  90.      * construct the rest of the Sturm sequence
  91.      */
  92.     for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
  93.  
  94.         /*
  95.          * reverse the sign and normalise
  96.          */
  97.         f = -fabs(sp->coef[sp->ord]);
  98.         for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  99.                 *fp /= f;
  100.     }
  101.  
  102.     sp->coef[0] = -sp->coef[0];    /* reverse the sign */
  103.  
  104.     return(sp - sseq);
  105. }
  106.  
  107. /*
  108.  * numroots
  109.  *
  110.  *    return the number of distinct real roots of the polynomial
  111.  * described in sseq.
  112.  */
  113. int
  114. numroots(np, sseq, atneg, atpos)
  115.         int        np;
  116.         poly    *sseq;
  117.         int        *atneg, *atpos;
  118. {
  119.         int        atposinf, atneginf;
  120.         poly    *s;
  121.         double    f, lf;
  122.  
  123.         atposinf = atneginf = 0;
  124.  
  125.  
  126.     /*
  127.      * changes at positve infinity
  128.      */
  129.     lf = sseq[0].coef[sseq[0].ord];
  130.  
  131.     for (s = sseq + 1; s <= sseq + np; s++) {
  132.             f = s->coef[s->ord];
  133.             if (lf == 0.0 || lf * f < 0)
  134.                 atposinf++;
  135.         lf = f;
  136.     }
  137.  
  138.     /*
  139.      * changes at negative infinity
  140.      */
  141.     if (sseq[0].ord & 1)
  142.             lf = -sseq[0].coef[sseq[0].ord];
  143.     else
  144.             lf = sseq[0].coef[sseq[0].ord];
  145.  
  146.     for (s = sseq + 1; s <= sseq + np; s++) {
  147.             if (s->ord & 1)
  148.                 f = -s->coef[s->ord];
  149.             else
  150.                 f = s->coef[s->ord];
  151.             if (lf == 0.0 || lf * f < 0)
  152.                 atneginf++;
  153.             lf = f;
  154.     }
  155.  
  156.     *atneg = atneginf;
  157.     *atpos = atposinf;
  158.  
  159.     return(atneginf - atposinf);
  160. }
  161.  
  162. /*
  163.  * numchanges
  164.  *
  165.  *    return the number of sign changes in the Sturm sequence in
  166.  * sseq at the value a.
  167.  */
  168. int
  169. numchanges(np, sseq, a)
  170.     int        np;
  171.     poly    *sseq;
  172.     double    a;
  173.  
  174. {
  175.     int        changes;
  176.     double    f, lf;
  177.     poly    *s;
  178.  
  179.     changes = 0;
  180.  
  181.     lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
  182.  
  183.     for (s = sseq + 1; s <= sseq + np; s++) {
  184.             f = evalpoly(s->ord, s->coef, a);
  185.             if (lf == 0.0 || lf * f < 0)
  186.                 changes++;
  187.             lf = f;
  188.     }
  189.  
  190.     return(changes);
  191. }
  192.  
  193. /*
  194.  * sbisect
  195.  *
  196.  *    uses a bisection based on the sturm sequence for the polynomial
  197.  * described in sseq to isolate intervals in which roots occur,
  198.  * the roots are returned in the roots array in order of magnitude.
  199.  */
  200. sbisect(np, sseq, min, max, atmin, atmax, roots)
  201.     int        np;
  202.     poly    *sseq;
  203.     double    min, max;
  204.     int        atmin, atmax;
  205.     double    *roots;
  206. {
  207.     double    mid;
  208.     int        n1, n2, its, atmid, nroot;
  209.  
  210.     if ((nroot = atmin - atmax) == 1) {
  211.  
  212.         /*
  213.          * first try a less expensive technique.
  214.          */
  215.         if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
  216.             return;
  217.  
  218.  
  219.         /*
  220.          * if we get here we have to evaluate the root the hard
  221.          * way by using the Sturm sequence.
  222.          */
  223.         for (its = 0; its < MAXIT; its++) {
  224.                 mid = (min + max) / 2;
  225.  
  226.                 atmid = numchanges(np, sseq, mid);
  227.  
  228.                 if (fabs(mid) > RELERROR) {
  229.                     if (fabs((max - min) / mid) < RELERROR) {
  230.                         roots[0] = mid;
  231.                         return;
  232.                     }
  233.                 } else if (fabs(max - min) < RELERROR) {
  234.                     roots[0] = mid;
  235.                     return;
  236.                 }
  237.  
  238.                 if ((atmin - atmid) == 0)
  239.                     min = mid;
  240.                 else
  241.                     max = mid;
  242.             }
  243.  
  244.         if (its == MAXIT) {
  245.                 fprintf(stderr, "sbisect: overflow min %f max %f\
  246.                     diff %e nroot %d n1 %d n2 %d\n",
  247.                     min, max, max - min, nroot, n1, n2);
  248.             roots[0] = mid;
  249.         }
  250.  
  251.         return;
  252.     }
  253.  
  254.     /*
  255.      * more than one root in the interval, we have to bisect...
  256.      */
  257.     for (its = 0; its < MAXIT; its++) {
  258.  
  259.             mid = (min + max) / 2;
  260.  
  261.             atmid = numchanges(np, sseq, mid);
  262.  
  263.             n1 = atmin - atmid;
  264.             n2 = atmid - atmax;
  265.  
  266.  
  267.             if (n1 != 0 && n2 != 0) {
  268.                 sbisect(np, sseq, min, mid, atmin, atmid, roots);
  269.                 sbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
  270.                 break;
  271.             }
  272.  
  273.             if (n1 == 0)
  274.                 min = mid;
  275.             else
  276.                 max = mid;
  277.     }
  278.  
  279.     if (its == MAXIT) {
  280.             fprintf(stderr, "sbisect: roots too close together\n");
  281.             fprintf(stderr, "sbisect: overflow min %f max %f diff %e\
  282.                 nroot %d n1 %d n2 %d\n",
  283.                 min, max, max - min, nroot, n1, n2);
  284.             for (n1 = atmax; n1 < atmin; n1++)
  285.             roots[n1 - atmax] = mid;
  286.     }
  287. }
  288.